//
// Created by chao.li on 2023/9/19.
//

#ifndef SRC_VEC_H
#define SRC_VEC_H

#include <cmath>
#include <iostream>

#pragma inline_depth( 255 )
#pragma inline_recursion( on )
#pragma auto_inline( on )

namespace vec_h
{
    struct data_xyz
    {
        double x, y, z;

        inline double &operator[](const int i)
        {
            return ((double *) this)[i];
        }

        inline double operator[](int i) const
        {
            return ((double *) this)[i];
        }
    };

    struct data_xy
    {
        double x, y;

        inline double &operator[](const int i)
        {
            return ((double *) this)[i];
        }

        inline double operator[](int i) const
        {
            return ((double *) this)[i];
        }
    };

    template<int dim, class data>
    struct vec_basic
            : public data
    {
    public:

        inline const double Evaluate(const int i) const
        { return (*((data *) this))[i]; }

        //////////////////// ASSIGMENT ////////////////////
        template<class ta>
        inline const vec_basic<dim, data> &operator=(const ta &Exp)
        {
            for (int I = 0; I < dim; I++)
                (*((data *) this))[I] = Exp.Evaluate(I);
            return *this;
        }

//    template< class ta >
//    inline vec_basic<dim, data> ( const ta& Exp )
//    {
//        for( int I=0; I<dim; I++ )
//            (*((data*)this))[I] = Exp.Evaluate( I );
//    }

        inline vec_basic<dim, data>()
        {
            for (int I = 0; I < dim; I++)
                (*((data *) this))[I] = 0;
        }

        template<class t_num0, class t_num1>
        inline vec_basic<dim, data>(t_num0 a, t_num1 b)
        {
            (*((data *) this))[0] = (double) a;
            (*((data *) this))[1] = (double) b;
        }

        template<class t_num0, class t_num1, class t_num2>
        inline vec_basic<dim, data>(t_num0 a, t_num1 b, t_num2 c)
        {
            (*((data *) this))[0] = (double) a;
            (*((data *) this))[1] = (double) b;
            (*((data *) this))[2] = (double) c;
        }

        template<class t_num0, class t_num1, class t_num2, class t_num3>
        inline vec_basic<dim, data>(t_num0 a, t_num1 b, t_num2 c, t_num3 d)
        {
            (*((data *) this))[0] = (double) a;
            (*((data *) this))[1] = (double) b;
            (*((data *) this))[2] = (double) c;
            (*((data *) this))[3] = (double) d;
        }

        template<class ta>
        inline vec_basic<dim, data>(const ta &Exp)
        {
            for (int I = 0; I < dim; I++)
                (*((data *) this))[I] = Exp.Evaluate(I);
        }

        inline double length() const
        {
            return sqrt(length_squared());
        }

        inline double length_squared() const
        {
            if (dim == 2) {
                return Evaluate(0) * Evaluate(0) + Evaluate(1) * Evaluate(1);
            } else if (dim == 3) {
                return Evaluate(0) * Evaluate(0) + Evaluate(1) * Evaluate(1) + Evaluate(2) * Evaluate(2);
            } else if (dim == 4) {
                return Evaluate(0) * Evaluate(0) + Evaluate(1) * Evaluate(1) + Evaluate(2) * Evaluate(2) +
                        Evaluate(3) * Evaluate(3);
            } else {
                return Evaluate(0) * Evaluate(0);
            }
        }

        inline bool near_zero() const
        {
            // Return true if the vector is close to zero in all dimensions.
            const auto s = 1e-8;
            if (dim == 2) {
                return fabs(Evaluate(0)) < s && fabs(Evaluate(1)) < s;
            } else if (dim == 3) {
                return fabs(Evaluate(0)) < s && fabs(Evaluate(1)) < s && fabs(Evaluate(2)) < s;
            } else if (dim == 4) {
                return fabs(Evaluate(0)) < s && fabs(Evaluate(1)) < s && fabs(Evaluate(2)) < s && fabs(Evaluate(3)) < s;
            } else {
                return fabs(Evaluate(0)) < s;
            }
        }
    };

//template<>
//struct vec_basic< 3, data_xyz>
//        : public data_xyz
//{
//    inline vec_basic< 3, data_xyz >( const double px, const double py, const double pz )
//    { x = px; y = py; z = pz; }
//};
//
    using vec3 = vec_basic<3, data_xyz>;
    using vec2 = vec_basic<2, data_xy>;

//////////////////// ARGUMENTS ////////////////////
    template<class ta_a>
    class vecarg
    {
        const ta_a &Argv;
    public:
        inline vecarg(const ta_a &A) : Argv(A)
        {}

        inline const double Evaluate(const int i) const
        { return Argv.Evaluate(i); }
    };

    template<>
    class vecarg<const double>
    {
        const double &Argv;
    public:
        inline vecarg(const double &A) : Argv(A)
        {}

        inline const double Evaluate(const int i) const
        { return Argv; }
    };

    template<>
    class vecarg<const float>
    {
        const float &Argv;
    public:
        inline vecarg(const float &A) : Argv(A)
        {}

        inline const float Evaluate(const int i) const
        { return Argv; }
    };

    template<>
    class vecarg<const int>
    {
        const int &Argv;
    public:
        inline vecarg(const int &A) : Argv(A)
        {}

        inline const double Evaluate(const int i) const
        { return (double) Argv; }
    };

//////////////////// EXPRESIONS ////////////////////
    template<class ta_a, class ta_b, class ta_eval>
    class vecexp_2
    {
        const vecarg<ta_a> Arg1;
        const vecarg<ta_b> Arg2;
    public:
        inline vecexp_2(const ta_a &A1, const ta_b &A2)
                : Arg1(A1), Arg2(A2)
        {}

        inline const double Evaluate(const int i) const
        { return ta_eval::Evaluate(i, Arg1, Arg2); }
    };

    template<class ta_a, class ta_eval>
    class vecexp_1
    {
        const vecarg<ta_a> Arg1;
    public:
        inline vecexp_1(const ta_a &A1) : Arg1(A1)
        {}

        inline const double Evaluate(const int i) const
        { return ta_eval::Evaluate(i, Arg1.Evaluate(i)); }
    };

//////////////////// SUM ////////////////////
    struct sum
    {
        template<class ta_a, class ta_b>
        inline static const double Evaluate(const int i, const ta_a &A, const ta_b &B)
        { return A.Evaluate(i) + B.Evaluate(i); }
    };

    struct sub
    {
        template<class ta_a, class ta_b>
        inline static const double Evaluate(const int i, const ta_a &A, const ta_b &B)
        { return A.Evaluate(i) - B.Evaluate(i); }
    };

    struct mul
    {
        template<class ta_a, class ta_b>
        inline static const double Evaluate(const int i, const ta_a &A, const ta_b &B)
        { return A.Evaluate(i) * B.Evaluate(i); }
    };

    struct divide
    {
        template<class ta_a, class ta_b>
        inline static const double Evaluate(const int i, const ta_a &A, const ta_b &B)
        { return A.Evaluate(i) / B.Evaluate(i); }
    };

    template<class ta_c1, class ta_c2>
    inline const vecexp_2<const ta_c1, const ta_c2, sum>
    operator+(const ta_c1 &Pa, const ta_c2 &Pb)
    { return vecexp_2<const ta_c1, const ta_c2, sum>(Pa, Pb); }

    template<class ta_c1, class ta_c2>
    inline const vecexp_2<const ta_c1, const ta_c2, sub>
    operator-(const ta_c1 &Pa, const ta_c2 &Pb)
    { return vecexp_2<const ta_c1, const ta_c2, sub>(Pa, Pb); }

    template<class ta_c1, class ta_c2>
    inline const vecexp_2<const ta_c1, const ta_c2, mul>
    operator*(const ta_c1 &Pa, const ta_c2 &Pb)
    { return vecexp_2<const ta_c1, const ta_c2, mul>(Pa, Pb); }

    template<class ta_c1, class ta_c2>
    inline const vecexp_2<const ta_c1, const ta_c2, divide>
    operator/(const ta_c1 &Pa, const ta_c2 &Pb)
    { return vecexp_2<const ta_c1, const ta_c2, divide>(Pa, Pb); }

    template<int dim, class data>
    inline vec_basic<dim, data> fabs(vec_basic<dim, data> v)
    {
        if (dim == 1) {
            return vec_basic<1, data>(fabs(v.x));
        } else if (dim == 2) {
            return vec_basic<2, data>(fabs(v.x), fabs(v.y));
        } else if (dim == 3) {
            return vec_basic<3, data>(fabs(v.x), fabs(v.y), fabs(v.z));
        } else if (dim == 4) {
            return vec_basic<4, data>(fabs(v.x), fabs(v.y), fabs(v.z), fabs(v.w));
        } else {
            return vec_basic<1, data>(fabs(v.x));
        }
    }

}
using vec3 = vec_h::vec3;
using vec2 = vec_h::vec2;

// vec3 Utility Functions
inline std::ostream &operator<<(std::ostream &out, const vec3 &v)
{
    return out << v[0] << ' ' << v[1] << ' ' << v[2];
}

inline std::ostream &operator<<(std::ostream &out, const vec2 &v)
{
    return out << v[0] << ' ' << v[1];
}

inline void test_vector()
{
    vec3 b(1, 2, 3), c(4, 5, 6);
    vec3 a = b + c;
    std::cout << a << std::endl;
}
inline double dot(const vec3 &u, const vec3 &v)
{
    return u[0] * v[0]
            + u[1] * v[1]
            + u[2] * v[2];
}

inline double dot(const vec2 &u, const vec2 &v)
{
    return u[0] * v[0]
            + u[1] * v[1];
}

inline vec3 cross(const vec3 &u, const vec3 &v)
{
    return vec3(u[1] * v[2] - u[2] * v[1],
                u[2] * v[0] - u[0] * v[2],
                u[0] * v[1] - u[1] * v[0]);
}

inline vec3 unit_vector(vec3 v)
{
    return v / v.length();
}

inline double frac(double v)
{
    return v > 0 ? v - floor(v) : v - ceil(v);
}

inline vec3 frac(vec3 v)
{
    return vec3(frac(v.x), frac(v.y), frac(v.z));
}

inline vec3 fabs(vec3 v)
{
    return vec3(fabs(v.x), fabs(v.y), fabs(v.z));
}

inline double clamp(double x, double min, double max)
{
    if (x < min) return min;
    if (x > max) return max;
    return x;
}

inline double saturate(double x)
{
    return clamp(x, 0, 1);
}

inline vec3 saturate(vec3 v)
{
    return vec3(saturate(v.x), saturate(v.y), saturate(v.z));
}

inline double lerp(double a, double b, double t)
{
    return a + (b - a) * t;
}

inline vec3 lerp(vec3 c1, vec3 c2, double t)
{
    return vec3(lerp(c1.x, c2.x, t),
                lerp(c1.y, c2.y, t),
                lerp(c1.z, c2.z, t));
}
#endif //SRC_VEC_H
